library(landscapemetrics)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggrepel)
library(readxl)
library(stars)
## Loading required package: abind
## Loading required package: sf
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
Восстанавливаем цветовую шкалу
typesdf = read_excel('data/types.xlsx')
types = typesdf$type
colormap = read_table2('data/colormap.clr',
col_names = c('N', 'R', 'G', 'B', 'A', 'type')) %>%
filter(type %in% types)
## Parsed with column specification:
## cols(
## N = col_double(),
## R = col_double(),
## G = col_double(),
## B = col_double(),
## A = col_double(),
## type = col_double()
## )
pal = rgb(colormap$R/255, colormap$G/255, colormap$B/255, colormap$A/255)
Грузим данные
read_lc = function(file) {
read_stars(file) %>%
select(type = 1) %>%
mutate(typefac = factor(type, levels = types),
urban = (type == 50),
crops = (type == 40)) %>%
st_warp(crs = 32662)
}
salekhard = read_lc('data/SALEKHARD.tif')
moscow = read_lc('data/MOSCOW.tif')
rostov = read_lc('data/ROSTOV.tif')
ufa = read_lc('data/UFA.tif')
petro = read_lc('data/PETRO.tif')
surgut = read_lc('data/SURGUT.tif')
ustug = read_lc('data/USTUG.tif')
grozny = read_lc('data/GROZNY.tif')
voronezh = read_lc('data/VORONEZH.tif')
Визуализируем:
plot_lc = function(lc) {
ggplot() +
geom_stars(data = lc['typefac']) +
scale_fill_manual(values = pal, breaks = types, labels = typesdf$class, drop = FALSE) +
coord_sf(crs = st_crs(lc)) +
xlab('LON') + ylab('LAT') +
theme_bw() + theme(legend.position = 'bottom', legend.title=element_blank()) +
guides(fill = guide_legend(ncol=4))
}
plot_lc(salekhard)
## Warning: Removed 4812 rows containing missing values (geom_raster).
plot_lc(moscow)
## Warning: Removed 4812 rows containing missing values (geom_raster).
plot_lc(rostov)
## Warning: Removed 4812 rows containing missing values (geom_raster).
plot_lc(ufa)
## Warning: Removed 4812 rows containing missing values (geom_raster).
plot_lc(petro)
## Warning: Removed 4812 rows containing missing values (geom_raster).
plot_lc(surgut)
## Warning: Removed 4812 rows containing missing values (geom_raster).
plot_lc(ustug)
## Warning: Removed 3208 rows containing missing values (geom_raster).
plot_lc(grozny)
## Warning: Removed 4812 rows containing missing values (geom_raster).
plot_lc(voronezh)
## Warning: Removed 4812 rows containing missing values (geom_raster).
Считаем метрики сложности:
lands = lst(САЛЕХАРД = salekhard, МОСКВА = moscow, `РОСТОВ-НА-ДОНУ` = rostov, УФА = ufa, ПЕТРОЗАВОДСК = petro,
ШАРЬЯ = ustug, ГРОЗНЫЙ = grozny, ВОРОНЕЖ = voronezh, СУРГУТ = surgut)
(metrics = imap(lands, function(land, name) {
tibble(name = name,
ent = lsm_l_ent(land) %>% pull(value),
condent = lsm_l_condent(land) %>% pull(value),
joinent = lsm_l_joinent(land) %>% pull(value),
mutinf = lsm_l_mutinf(land) %>% pull(value),
relmutinf = mutinf / ent,
urban = 100 * sum(land$urban, na.rm = TRUE) / length(land$urban),
crops = 100 * sum(land$crops, na.rm = TRUE) / length(land$crops))
}) %>% bind_rows())
## # A tibble: 9 x 8
## name ent condent joinent mutinf relmutinf urban crops
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 САЛЕХАРД 2.52 1.13 3.65 1.39 0.552 0.115 0
## 2 МОСКВА 2.95 1.38 4.33 1.57 0.532 17.4 17.3
## 3 РОСТОВ-НА-ДОНУ 1.49 0.561 2.05 0.927 0.623 3.20 72.4
## 4 УФА 2.16 0.977 3.13 1.18 0.547 3.48 45.9
## 5 ПЕТРОЗАВОДСК 2.35 0.945 3.30 1.41 0.599 0.196 0.985
## 6 ШАРЬЯ 2.03 0.853 2.88 1.18 0.579 0.245 4.47
## 7 ГРОЗНЫЙ 2.16 0.780 2.94 1.38 0.640 2.75 22.6
## 8 ВОРОНЕЖ 1.86 0.759 2.62 1.10 0.592 3.09 65.5
## 9 СУРГУТ 2.90 1.27 4.17 1.63 0.561 0.343 0.0543
ggplot(metrics, mapping = aes(ent, relmutinf, label = name)) +
geom_point()
Читаем данные:
scales = c(200, 500, 1000)
regions = c('МОСКВА', 'ПЕТРОЗАВОДСК', 'САЛЕХАРД', 'УФА', 'РОСТОВ-ДОН')
lindata = lapply(scales, function(scale) {
read_delim(paste0('tables/resLines', scale, '.txt'), delim = ';', locale = locale(decimal_mark = ".")) %>%
select(-ncol(.)) %>%
mutate(Scale = scale,
Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(as.character(Scale))))
}) %>%
bind_rows() %>%
filter(!stringr::str_detect(Name, 'clip')) %>%
mutate(Layer = stringr::str_sub(Name, 1, 3),
Dim = stringr::str_sub(Name, 4, 6)) %>%
filter(!(Layer %in% c('rlh', 'rlf', 'veg'))) %>%
mutate_all(~replace(., is.nan(.), 0))
## Warning: Missing column names filled in: 'X17' [17]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## PointsNumber = col_double(),
## BendNumber = col_double(),
## AveAngle = col_double(),
## TotalArea = col_double(),
## Length = col_double(),
## AverageLength = col_double(),
## AverageBaseLineLength = col_double(),
## AverageHeight = col_double(),
## AverageWidth = col_double(),
## AverageArea = col_double(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X17 = col_logical()
## )
## Warning: Missing column names filled in: 'X17' [17]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## PointsNumber = col_double(),
## BendNumber = col_double(),
## AveAngle = col_double(),
## TotalArea = col_double(),
## Length = col_double(),
## AverageLength = col_double(),
## AverageBaseLineLength = col_double(),
## AverageHeight = col_double(),
## AverageWidth = col_double(),
## AverageArea = col_double(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X17 = col_logical()
## )
## Warning: Missing column names filled in: 'X17' [17]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## PointsNumber = col_double(),
## BendNumber = col_double(),
## AveAngle = col_double(),
## TotalArea = col_double(),
## Length = col_double(),
## AverageLength = col_double(),
## AverageBaseLineLength = col_double(),
## AverageHeight = col_double(),
## AverageWidth = col_double(),
## AverageArea = col_double(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X17 = col_logical()
## )
pntdata = lapply(scales, function(scale) {
read_delim(paste0('tables/resPoints', scale, '.txt'), delim = ';', locale = locale(decimal_mark = ".")) %>%
select(-ncol(.)) %>%
mutate(Scale = scale,
Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(as.character(Scale))))
}) %>%
bind_rows() %>%
mutate(Layer = stringr::str_sub(Name, 1, 3),
Dim = stringr::str_sub(Name, 4, 6)) %>%
filter(!(Layer %in% c('rlh', 'rlf', 'veg')))
## Warning: Missing column names filled in: 'X7' [7]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X7 = col_logical()
## )
## Warning: Missing column names filled in: 'X7' [7]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X7 = col_logical()
## )
## Warning: Missing column names filled in: 'X7' [7]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X7 = col_logical()
## )
poldata = lapply(scales, function(scale) {
read_delim(paste0('tables/resPolygons', scale, '.txt'), delim = ';', locale = locale(decimal_mark = ".")) %>%
select(-ncol(.)) %>%
mutate(Scale = scale,
Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(as.character(Scale))))
}) %>%
bind_rows() %>%
mutate(Layer = stringr::str_sub(Name, 1, 3),
Dim = stringr::str_sub(Name, 4, 6)) %>%
filter(!(Layer %in% c('rlh', 'rlf', 'veg')))
## Warning: Missing column names filled in: 'X12' [12]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## PointsNumber = col_double(),
## TotalArea = col_double(),
## TotalPerimetr = col_double(),
## AveArea = col_double(),
## AvePerimetr = col_double(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X12 = col_logical()
## )
## Warning: Missing column names filled in: 'X12' [12]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## PointsNumber = col_double(),
## TotalArea = col_double(),
## TotalPerimetr = col_double(),
## AveArea = col_double(),
## AvePerimetr = col_double(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X12 = col_logical()
## )
## Warning: Missing column names filled in: 'X12' [12]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Name = col_character(),
## PointsNumber = col_double(),
## TotalArea = col_double(),
## TotalPerimetr = col_double(),
## AveArea = col_double(),
## AvePerimetr = col_double(),
## AtrbCount = col_double(),
## UniqueValues = col_double(),
## AveUniqueValues = col_double(),
## NumberOfObj = col_double(),
## X12 = col_logical()
## )
nuniq = length(unique(poldata$Region))
inters = read_delim('tables/topLines.txt', delim = ';', locale = locale(decimal_mark = ".")) %>%
select(-ncol(.)) %>%
mutate(Scale = as.character(rep(scales, each = n()/3)),
Region = stringr::str_sub(Region, 1, nchar(Region) - nchar(Scale)))
## Warning: Missing column names filled in: 'X3' [3]
## Parsed with column specification:
## cols(
## Region = col_character(),
## Intersections = col_double(),
## X3 = col_logical()
## )
Подсчитаем суммарное количество социально-экономических точек:
linpts = lindata %>%
group_by(Region, Scale) %>%
summarise(linnpts = sum(PointsNumber))
polpts = poldata %>%
group_by(Region, Scale) %>%
summarise(polnpts = sum(PointsNumber),
area = 1e-6 * sum(TotalArea * (Layer == 'adm')))
pntpts = pntdata %>%
group_by(Region, Scale) %>%
summarise(pntnpts = sum(NumberOfObj))
stats = bind_cols(linpts, polpts, pntpts) %>%
select(name = Region, Scale, pntnpts, linnpts, polnpts, area) %>%
mutate(total = pntnpts + linnpts + polnpts,
density = total / area) %>%
ungroup() %>%
mutate(name = ifelse(name == 'РОСТОВ-ДОН', 'РОСТОВ-НА-ДОНУ',
ifelse(name == 'УСТЮГ', 'ШАРЬЯ', name)))
Соединяем
compare = left_join(stats, metrics)
## Joining, by = "name"
Корреляционный и регрессионный анализ
tab = compare %>%
filter(Scale == 1000)
ggplot(tab, mapping = aes(urban, density, label = name)) +
geom_point() +
geom_smooth(method = 'lm') +
geom_text_repel(
segment.size = 0.3,
segment.color = "grey50",
direction = "y",
hjust = 0
) +
scale_x_log10()
ggplot(tab, mapping = aes(crops, density, label = name)) +
geom_point() +
geom_smooth(method = 'lm') +
geom_text_repel(
segment.size = 0.3,
segment.color = "grey50",
direction = "y",
hjust = 0
) +
scale_x_log10()
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
ggplot(tab, mapping = aes(crops + urban, density, label = name)) +
geom_point() +
geom_smooth(method = 'lm') +
geom_text_repel(
segment.size = 0.3,
segment.color = "grey50",
direction = "y",
hjust = 0
) +
scale_x_log10()
cor.test(tab$density, log(tab$urban + 1))
##
## Pearson's product-moment correlation
##
## data: tab$density and log(tab$urban + 1)
## t = 4.9167, df = 7, p-value = 0.00172
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5214060 0.9746945
## sample estimates:
## cor
## 0.8805971
cor.test(tab$density, log(tab$crops + 1))
##
## Pearson's product-moment correlation
##
## data: tab$density and log(tab$crops + 1)
## t = 3.0497, df = 7, p-value = 0.01859
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1830774 0.9452803
## sample estimates:
## cor
## 0.7553583
cor.test(tab$density, log(tab$crops + 1) + log(tab$urban + 1))
##
## Pearson's product-moment correlation
##
## data: tab$density and log(tab$crops + 1) + log(tab$urban + 1)
## t = 4.3883, df = 7, p-value = 0.003202
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4458400 0.9692525
## sample estimates:
## cor
## 0.8563911
lm(density ~ log(urban + 1), tab)
##
## Call:
## lm(formula = density ~ log(urban + 1), data = tab)
##
## Coefficients:
## (Intercept) log(urban + 1)
## 0.0646 0.1251
lm(density ~ log(crops + 1), tab)
##
## Call:
## lm(formula = density ~ log(crops + 1), data = tab)
##
## Coefficients:
## (Intercept) log(crops + 1)
## 0.06255 0.05715
lm(density ~ log(crops + urban + 1), tab)
##
## Call:
## lm(formula = density ~ log(crops + urban + 1), data = tab)
##
## Coefficients:
## (Intercept) log(crops + urban + 1)
## 0.04241 0.06152
Зависимость плотности от сложности ландшафта:
ggplot(tab, mapping = aes(joinent, density, label = name)) +
geom_point() +
geom_text_repel(
segment.size = 0.3,
segment.color = "grey50",
direction = "y",
hjust = 0
)